# PACOTES -----------------------------------------------------------------
pacman::p_load(pacman, ggplot2, tidyverse, magick)
# FUNÇÃO ------------------------------------------------------------------
# Esta função faz os cálculos de transiente hidráulico em uma adutora assentada em z = 0 m
# somente com duas condições de contorno: uma imediatamente à montante caracterizada por um reservatório
# de nível fixo na posição x = 0 m; e uma válvula localizada no final do comprimento da adutora com fechamento
# rápido com tc = 2l/a.
fun.transiente <- function(na.reservatorio, # cota do nível d'água (m)
l.adutora, # comprimento da adutora (m)
d.adutora, # diâmetro da adutora (m)
f.atrito, # fator de atrito (fórmula universal)
celeridade, # celeriadade da onda de pressão (m/s)
vel.permanente, # velocidade em regime permanente (m/s)
tempo.simulacao, # tempo de simulação (s)
passo.x # discretização horizontal (m)
){
# Início contador de tempo
tempo.inicio <- Sys.time()
# Pacotes
if (!require("pacman")) install.packages("pacman")
p_load(pacman, tidyverse, beepr)
# Renomear parâmetros
h.r <- na.reservatorio
l <- l.adutora
d <- d.adutora
f <- f.atrito
a <- celeridade
v0 <- vel.permanente
t.s <- tempo.simulacao
dx <- passo.x
g <- 9.81
# Discretização
n.dx <- l/dx # número de passos na horizontal
dt <- dx/a # discretização temporal
n.dt <- t.s/dt # número de passos no tempo
# Inicialização das variáveis Vazão e Carga
q <- matrix(0, ncol = n.dx, nrow = n.dt) %>% data.frame # vazão (m³/s)
h <- matrix(0, ncol = n.dx, nrow = n.dt) %>% data.frame # carga de pressão (m)
# Perda de carga unitária em cada passo dx
dh <- f*l/d*v0^2/(2*g)/n.dx
# Constantes
area <- pi*(d/2)^2 # área (m²)
b <- a/(g*area) # impedância característica
tc <- 2*l/a # tempo de fechamento da válvula (s)
r <- f*dx/(2*g*d*area^2) # coeficiente de resistência
# Regime permanente (t = 0) e condições de contorno
q[1,] <- v0*area # primeira linha (t = 0) calculada a partir da vazão em regime permanente
h[,1] <- h.r # primeira coluna (x = 0) é o nível d'água no reservatório
for(i in 2:n.dx){ # propagar a perda de carga ao longo da adutora
h[1,i] <- h[1,i-1] - dh
}
# Regime transitório
for(t in 2:n.dt){
# Contador
# cat(paste("\npasso de tempo atual:", (t-1)*dt, "segundos.\n"))
# Fechamento da válvula
tempo <- (t-1)*dt
if(tempo <= 3){
tau <- (1 - tempo/tc)^3
}
if(tempo > 3){
tau <- 0
}
cv <- (q[1,n.dx]*tau)^2/(2*h[1,n.dx])
# Iterações ao longo da adutora
for(i in 1:n.dx){
# Condições de montante (reservatório)
if(i == 1){
cb <- h[t-1, i+1] - b*q[t-1, i+1]
bb <- b + r*abs(q[t-1, i+1])
q[t,i] <- (h[t,i] - cb)/bb
}
if(i > 1 & i < n.dx){
ca <- h[t-1, i-1] + b*q[t-1, i-1]
ba <- b + r*abs(q[t-1, i-1])
cb <- h[t-1, i+1] - b*q[t-1, i+1]
bb <- b + r*abs(q[t-1, i+1])
h[t,i] <- (bb*ca + ba*cb)/(ba + bb) # carga no ponto P em xi
q[t,i] <- (ca - cb)/(ba + bb) # vazão no ponto P em xi
}
if(i == n.dx){
ca <- h[t-1, i-1] + b*q[t-1, i-1]
ba <- b + r*abs(q[t-1, i-1])
# Vazão e carga na válvula
q[t,i] <- -ba*cv + sqrt((ba*cv)^2 + 2*cv*ca)
h[t,i] <- ca - ba*q[t,i]
}
}
}
# Adicionar coluna de 'tempo'
t <- seq(from = 0, to = t.simulacao, by = dt) %>% head(-1) # remove a última linha a mais criada
h <- cbind(t, h)
q <- cbind(t, q)
resultados <- list(Carga = h, Vazao = q)
# Alerta sonoro
beepr::beep(sound = 10)
tempo.decorrido <- difftime(Sys.time(), tempo.inicio, units = "secs")
cat("\nSimulação de", t.s, "segundos concluída em:", tempo.decorrido,"segundos\n")
return(resultados)
}
# DADOS DE ENTRADA --------------------------------------------------------
na.reservatorio <- 100 # [m] nível d'água do reservatório de nível fixo
l.adutora <- 1500 # [m] extensão da adutora
d.adutora <- 1 # [m] diâmetro interno da adutora
f.atrito <- 0.02 # fator de atrito para condições de escoamento permanente e transitório
na.adutora <- 0 # [m] cota horizontal da adutora
v0 <- 1 # [m/s] velocidade de escoamento em regime permanente
a <- 1000 # [m/s] celeridade da onda de pressão
g <- 9.81 # [m/s²] aceleração da gravidade
t.simulacao <- 60 # [s] tempo de simulação
# DISCRETIZAÇÃO ---------------------------------------------------------------
# Espacial [m]
dx <- 10
# SIMULAÇÃO ---------------------------------------------------------------
# Rodar função p/ simular transiente
resultado <- fun.transiente(na.reservatorio = 100,
l.adutora = 1500,
d.adutora = 1,
f.atrito = 0.02,
celeridade = 1000,
vel.permanente = 1,
tempo.simulacao = 60,
passo.x = 10)
##
## Simulação de 60 segundos concluída em: 141.9613 segundos
# VISUALIZAÇÃO DOS RESULTADOS ---------------------------------------------
# Plotar a válvula (coluna 'X150') durante os primeiros 60 segundos de simulação
plot.h.valvula <- ({
ggplot() +
geom_line(data = resultado$Carga, aes(x = t, y = X150), color = "steelblue", linewidth = 0.4) +
geom_hline(yintercept = resultado$Carga[1,151],
color = "red", linetype = "dashed", linewidth = 0.4) +
labs(x = "Tempo (s)", y = "Carga de Pressão (m)") +
theme_minimal() +
theme(panel.backgroundd = NULL,
plot.background = element_rect(color = "white", fill = "white"),
panel.border = element_rect(fill = NA),
text = element_text(size = 10))
})
ggsave(filename = "Plotagens/Grafico A_Carga de Pressão na Válvula.png",
plot = plot.h.valvula,
height = 10, width = 16, units = "cm", dpi = 100)
plot.h.valvula
plot.q.valvula <- ({
ggplot() +
geom_line(data = resultado$Vazao, aes(x = t, y = X150), color = "steelblue", linewidth = 0.4) +
geom_hline(yintercept = resultado$Vazao[1,151],
color = "red", linetype = "dashed", linewidth = 0.4) +
labs(x = "Tempo (s)", y = "Vazão (m³/s)") +
theme_minimal() +
theme(panel.backgroundd = NULL,
plot.background = element_rect(color = "white", fill = "white"),
panel.border = element_rect(fill = NA),
text = element_text(size = 10))
})
ggsave(filename = "Plotagens/Grafico A_Vazão na Válvula.png",
plot = plot.q.valvula,
height = 10, width = 16, units = "cm", dpi = 100)
plot.q.valvula
# Plotar Carga de Pressão e Vazão em X = 750 m
plot.h.750 <- ({
ggplot() +
geom_line(data = resultado$Carga, aes(x = t, y = X75), color = "steelblue", linewidth = 0.4) +
geom_hline(yintercept = resultado$Carga[1,76],
color = "red", linetype = "dashed", linewidth = 0.4) +
labs(x = "Tempo (s)", y = "Carga de Pressão (m)") +
theme_minimal() +
theme(panel.backgroundd = NULL,
plot.background = element_rect(color = "white", fill = "white"),
panel.border = element_rect(fill = NA),
text = element_text(size = 10))
})
ggsave(filename = "Plotagens/Grafico A_Carga de Pressão na Válvula.png",
plot = plot.h.750,
height = 10, width = 16, units = "cm", dpi = 100)
plot.h.750
plot.q.750 <- ({
ggplot() +
geom_line(data = resultado$Vazao, aes(x = t, y = X75), color = "steelblue", linewidth = 0.4) +
geom_hline(yintercept = resultado$Vazao[1,76],
color = "red", linetype = "dashed", linewidth = 0.4) +
labs(x = "Tempo (s)", y = "Vazão (m³/s)") +
theme_minimal() +
theme(panel.backgroundd = NULL,
plot.background = element_rect(color = "white", fill = "white"),
panel.border = element_rect(fill = NA),
text = element_text(size = 10))
})
ggsave(filename = "Plotagens/Grafico A_Carga de Pressão na Válvula.png",
plot = plot.q.750,
height = 10, width = 16, units = "cm", dpi = 100)
plot.q.750
# VISUALIZAÇÃO ANIMADA ----------------------------------------------------
# Transformar dado wide em long
h.long <- ({
resultado$Carga %>%
pivot_longer(cols = -t,
names_to = "comprimento",
values_to = "carga") %>%
mutate(comprimento = as.numeric(gsub("X", "", comprimento))*dx) %>%
rename("tempo" = "t")
})
q.long <- ({
resultado$Vazao %>%
pivot_longer(cols = -t,
names_to = "comprimento",
values_to = "vazao") %>%
mutate(comprimento = as.numeric(gsub("X", "", comprimento))*dx) %>%
rename("tempo" = "t")
})
# Definir diretório temporário
dir.out.h <- file.path(tempdir(), "plot.carga")
dir.create(dir.out.h, recursive = TRUE)
h.min <- min(h.long$carga)
h.max <- max(h.long$carga)
# Loop p/ gerar frames
for(i in seq(from = 0, max(h.long$tempo), 0.1)){
p <- h.long %>%
# filter(tempo == i) %>% # remover o filtro de tempo
filter(abs(tempo - i) < 1e-6) %>%
ggplot(aes(x = comprimento, y = carga)) +
geom_line() +
geom_hline(yintercept = 100, color = "red", alpha = 0.75,
linewidth = 0.3, linetype = "dashed") +
scale_y_continuous(limits = c(h.min, h.max)) +
labs(x = "Comprimento da Adutora [m]", y = "Carga de Pressão [m]",
title = paste("Tempo:", i, "segundos")) +
theme_minimal() +
theme(panel.background = NULL,
plot.background = element_rect(color = "white", fill = "white"),
panel.border = element_rect(fill = NA),
text = element_text(size = 10))
fp <- file.path(dir.out.h, paste0("h_t_", i*10,".jpg"))
ggsave(plot = p,
filename = fp,
height = 12, width = 16, units = 'cm',
dpi = 75)
}
# Criar e baixar GIF
imgs.h <- list.files(dir.out.h, full.names = TRUE)
imgs.h <- imgs.h[order(as.numeric(gsub("[^0-9]", "", basename(imgs.h))))] # extrair frames na ordem
imgs.list.h <- lapply(imgs.h, image_read) # importar frames de volta p/ R
img.join.h <- image_join(imgs.list.h) # unir frames
img.animate.h <- image_animate(img.join.h, fps = 10) # gerar animação
image_write(image = img.animate.h, # criar arquivo GIF
path = "carga.gif")
img.animate.h
